function fig3_22

%  plots exact solution elliptic equation

ep=0.05; b=0;  a=1;

% set parameters
%nx=40; nn=40;
nx=40; nn=40;
x=linspace(-1,1,nx);
y=x;

alpha=1;  beta=1;
chi=0.5*sqrt(alpha^2 + beta^2 - 4*ep)/ep;

pi2=2*pi;
b0 = (b-a)*ifc(0,pi2,0,ep)/pi2;
bess0 = besseli(0,chi);
for in=1:nn
    ap(in)=(b-a)*ifs(0,pi2,in,ep)/pi;
    bp(in)=(b-a)*ifc(0,pi2,in,ep)/pi;
    bess(in)=besseli(in,chi);
end

U=zeros(nx,nx);
for ix=1:nx
    for iy=1:nx
        r=sqrt(x(ix)^2+y(iy)^2);
        if r < 1
            phi=atan2(y(iy),x(ix));
            rp=chi*r;
            s=b0*besseli(0,rp)/bess0;
            for in=1:nn
                t1=besseli(in,rp)*(ap(in)*sin(in*phi)+bp(in)*cos(in*phi))/bess(in);
                s=s+t1;
            end
            U(ix,iy)=a+s*exp(-0.5*(x(ix)+y(iy))/ep);
        end
    end
end

clf
% get(gcf)
%set(gcf,'Position', [7 1004 573 333])
set(gcf,'Position', [7 979 601 358])
co = [0 0 1;
    0 0.5 0;
    1 0 0;
    0 0.75 0.75;
    0.75 0 0.75;
    0.75 0.75 0;
    0.25 0.25 0.25];
set(groot,'defaultAxesColorOrder',co)
%subaxis(1,1,1,1,'MT',0.01,'MB',0.18,'MR',-0.02,'ML',0.05,'P',0.04)

UU=U+3;
hold on
grid on
box on
%axes('position',[.13  .43  .8  .53])

surf(x,y,UU','FaceColor','interp')
%surfc(x,y,U')
view(38, 22);
%colormap gray
colormap parula

% get(gca)
set(gca,'ztick',[0 1 2 3]);
set(gca,'ZTickLabel',{'-3';'-2';'-1';'0'})

set(gca,'xtick',[-1 0 1]);
set(gca,'ytick',[-1 0 1]);

% countour plots

nx=200; nn=nx;
x=linspace(-1,1,nx);
y=x;

alpha=1;  beta=1;
chi=0.5*sqrt(alpha^2 + beta^2 - 4*ep)/ep;

pi2=2*pi;
b0 = (b-a)*ifc(0,pi2,0,ep)/pi2;
bess0 = besseli(0,chi);
for in=1:nn
    ap(in)=(b-a)*ifs(0,pi2,in,ep)/pi;
    bp(in)=(b-a)*ifc(0,pi2,in,ep)/pi;
    bess(in)=besseli(in,chi);
end

U=zeros(nx,nx);
for ix=1:nx
    for iy=1:nx
        r=sqrt(x(ix)^2+y(iy)^2);
        if r < 1
            phi=atan2(y(iy),x(ix));
            rp=chi*r;
            s=b0*besseli(0,rp)/bess0;
            for in=1:nn
                t1=besseli(in,rp)*(ap(in)*sin(in*phi)+bp(in)*cos(in*phi))/bess(in);
                s=s+t1;
            end
            U(ix,iy)=a+s*exp(-0.5*(x(ix)+y(iy))/ep);
        end
    end
end

UU=U+3;
contour(x,y,UU',9,'k','LineWidth',0.4)
xlabel('$x$\,-\,axis','Interpreter','latex','FontSize',14)
ylabel('$y$\,-\,axis','Interpreter','latex','FontSize',14)
%zlabel('Solution','FontSize',14,'FontWeight','bold')
% Set the fontsize to 12 for the plot  '

ax = gca;
ax.XAxis.TickLabelInterpreter = 'latex';
ax.YAxis.TickLabelInterpreter = 'latex';
set(gca,'FontSize',13)

hold off

%exportgraphics(gcf,'/Users/mark/Desktop/fig216.pdf')

function f=fs(phi,n,ep)
f = exp(0.5*(cos(phi)+sin(phi))/ep)*sin(n*phi);

function f=fc(phi,n,ep)
f = exp(0.5*(cos(phi)+sin(phi))/ep)*cos(n*phi);

function s=ifs(a,b,n,ep)
%  simpson's method
N=max(100,100*n);
xd=linspace(a,b,N+1);
h=xd(2)-xd(1);
s=fs(xd(1),n,ep);
for j=2:2:N
    f2=fs(xd(j+1),n,ep);
    s=s+4*fs(xd(j),n,ep)+2*f2;
end
s=h*(s-f2)/3;

function s=ifc(a,b,n,ep)
%  simpson's method
N=max(100,100*n);
xd=linspace(a,b,N+1);
h=xd(2)-xd(1);
s=fc(xd(1),n,ep);
for j=2:2:N
    f2=fc(xd(j+1),n,ep);
    s=s+4*fc(xd(j),n,ep)+2*f2;
end
s=h*(s-f2)/3;















